{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scalable GP Regression (w/ KISS-GP)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "If the function you are modeling has additive structure across its dimensions, then SKI can be one of the most efficient methods for your problem.\n",
    "\n",
    "Here, we model the kernel as a sum of kernels that each act on one dimension. Additive SKI (or KISS-GP) can work very well out-of-the-box on larger datasets (100,000+ data points) with many dimensions. This is a strong assumption though, and may not apply to your problem.\n",
    "\n",
    "This is the same as [the KISSGP Kronecker notebook](./KISSGP_Kronecker_Regression.ipynb), but applied to more dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Inline plotting\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set up train data\n",
    "\n",
    "Here we're learning a simple sin function - but in 2 dimensions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We store the data as a 10k 1D vector\n",
    "# It actually represents [0,1]x[0,1] in cartesian coordinates\n",
    "n = 30\n",
    "train_x = torch.zeros(pow(n, 2), 2)\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        # Each coordinate varies from 0 to 1 in n=100 steps\n",
    "        train_x[i * n + j][0] = float(i) / (n-1)\n",
    "        train_x[i * n + j][1] = float(j) / (n-1)\n",
    "\n",
    "train_x = train_x.cuda()\n",
    "train_y = torch.sin(train_x[:, 0]) + torch.cos(train_x[:, 1]) * (2 * math.pi) + torch.randn_like(train_x[:, 0]).mul(0.01)\n",
    "train_y = train_y / 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model\n",
    "\n",
    "As with the Kronecker example case, applying SKI to a multidimensional kernel is as simple as wrapping that kernel with a `GridInterpolationKernel`. You'll want to be sure to set `num_dims` though!\n",
    "\n",
    "To use an additive decomposition of the kernel, wrap it in an `AdditiveStructureKernel`.\n",
    "\n",
    "SKI has only one hyperparameter that you need to worry about: the grid size. For 1D functions, a good starting place is to use as many grid points as training points. (Don't worry - the grid points are really cheap to use, especially with an additive function!).\n",
    "\n",
    "If you want, you can also explicitly determine the grid bounds of the SKI approximation using the `grid_bounds` argument. However, it's easier if you don't use this argument - then GPyTorch automatically chooses the best bounds for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPRegressionModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)\n",
    "        \n",
    "        # SKI requires a grid size hyperparameter. This util can help with that\n",
    "        # We're setting Kronecker structure to False because we're using an additive structure decomposition\n",
    "        grid_size = gpytorch.utils.grid.choose_grid_size(train_x, kronecker_structure=False)\n",
    "        \n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.AdditiveStructureKernel(\n",
    "            gpytorch.kernels.GridInterpolationKernel(\n",
    "                gpytorch.kernels.ScaleKernel(\n",
    "                    gpytorch.kernels.RBFKernel(),\n",
    "                ), grid_size=128, num_dims=1\n",
    "            ), num_dims=2\n",
    "        )\n",
    "\n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)\n",
    "\n",
    "    \n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood()\n",
    "model = GPRegressionModel(train_x, train_y, likelihood).cuda()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model hyperparameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/gpleiss/workspace/gpytorch/gpytorch/kernels/kernel.py:278: UserWarning: You are using a version of PyTorch where torch.pdist does not support batch matrices. Falling back on manual distance computation. Updating PyTorch to the latest pytorch-nightly build will offer significant memory savings during kernel computation.\n",
      "  warnings.warn('You are using a version of PyTorch where torch.pdist does not support batch '\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter 1/20 - Loss: 0.752\n",
      "Iter 2/20 - Loss: 0.715\n",
      "Iter 3/20 - Loss: 0.676\n",
      "Iter 4/20 - Loss: 0.638\n",
      "Iter 5/20 - Loss: 0.598\n",
      "Iter 6/20 - Loss: 0.559\n",
      "Iter 7/20 - Loss: 0.517\n",
      "Iter 8/20 - Loss: 0.477\n",
      "Iter 9/20 - Loss: 0.434\n",
      "Iter 10/20 - Loss: 0.393\n",
      "Iter 11/20 - Loss: 0.351\n",
      "Iter 12/20 - Loss: 0.307\n",
      "Iter 13/20 - Loss: 0.262\n",
      "Iter 14/20 - Loss: 0.218\n",
      "Iter 15/20 - Loss: 0.172\n",
      "Iter 16/20 - Loss: 0.127\n",
      "Iter 17/20 - Loss: 0.081\n",
      "Iter 18/20 - Loss: 0.033\n",
      "Iter 19/20 - Loss: -0.015\n",
      "Iter 20/20 - Loss: -0.061\n",
      "CPU times: user 1.41 s, sys: 140 ms, total: 1.55 s\n",
      "Wall time: 1.55 s\n"
     ]
    }
   ],
   "source": [
    "# Optimize the model\n",
    "model.train()\n",
    "likelihood.train()\n",
    "\n",
    "# Use the adam optimizer\n",
    "optimizer = torch.optim.Adam([\n",
    "    {'params': model.parameters()},  # Includes GaussianLikelihood parameters\n",
    "], lr=0.1)\n",
    "\n",
    "# \"Loss\" for GPs - the marginal log likelihood\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "# Sometimes we get better performance on the GPU when we don't use Toeplitz math\n",
    "# for SKI. This flag controls that\n",
    "def train(num_iter):\n",
    "    with gpytorch.settings.use_toeplitz(False):\n",
    "        for i in range(num_iter):\n",
    "            optimizer.zero_grad()\n",
    "            output = model(train_x)\n",
    "            loss = -mll(output, train_y)\n",
    "            loss.backward()\n",
    "            print('Iter %d/%d - Loss: %.3f' % (i + 1, num_iter, loss.item()))\n",
    "            optimizer.step()\n",
    "            \n",
    "%time train(num_iter=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO8AAADSCAYAAACilJfKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXm0bEd1n79f3/eEmIQkHmBNIDCyA2ZGIBLwQkwGxJgsIAgCgjCEFWJgGcJshjAYbAfjLIxlBbAQGDFjZALIipAgGDMIA2JQGA3ogYKsiUGAdO/tnT/Oue/uqu6qrr5Td7+3v7V63XNOjd23q2vXrr1ry8wIgmDxGMy6A0EQbIwYvEGwoMTgDYIFJQZvECwoMXiDYEGJwRsEC8pCDl5Jx0oySbv6+49JOmUH2n2FpHcuWt0T2v0jSc9pyHeqpD/sr0+UtHeD7X1f0v376xdLekt/nfxPtxNJ50t6an/9cEnv3u42t4NtG7z9P+lXkn4h6SeS/lrSDbajLTN7sJm9vbFP99/q9iUdJWlF0m+OSfuQpD/d6ja3Akk3AZ4I/FV/XxyUZvYMM3vVVrZvZq81s6duZZ0b6MNZwO0k3WGW/dgI2z3zPszMbgDcBbgb8NI8gzoWUgJYw8x+BJwLPME/l3Q4cBIw8YdlRjwJ+KiZ/WrWHZkxZwJPn3UnpmVHBk3/5f4YcDvYJ7a8RtI/AL8EbiXpRpLeKukSST+S9GpJS33+JUl/KukySd8DHuLr92JQf/80SRdJ+rmkb0i6i6R3ADcH/q6XBp7f572HpM9IukrSVySd6Oq5paRP9vWcA+ypvM23kw1e4LHA183sq319fy7pYkk/k/RFSb87rqJxM2Ambg4kvVDSdyVdLum9/Q8Fkg6W9M7++VWSviDpZoU+Pxj4ZOU9+fZPl/TqQtqz+s/56P7+oZK+3Lf/mdKsVlgqPF7SD/v/9Utc3utIeqOkH/evN0q6jkt/mqTvSLpC0lmSjnRpD5D0fyX9VNKbAGVtnk/2nVoEdmTwSjqGbgb6knv8BLpfuxsCP6D78q8AtwbuDPwesDYgnwY8tH9+PPCoSluPBl5BJw4eAjwcuNzMngD8kF4aMLM/lnQU8L+AVwOHA88DPtCLkwDvAr5IN2hfBdTW1R8C9ki6V/Yez3D3XwDu1Lf1LuB9kg6u1FniWcAjgXsDRwJXAn/Rp50C3Ag4Brgx8AygNLPeHvjmBtrfR78OfhJwbzPbK+kuwNuA/9S3/1fAWX6gTeBewG8D9wNeJuk2/fOXAPeg+/zuCNydXpKTdF/gj4DHAEfQfZ/e3aftAT7Q590DfBe4Z9bmRcCxkg6Z5r3PHDPblhfwfeAXwFV0H+abgev2aecD/83lvRlwzVp6/+xk4Lz++hPAM1za7wEG7HL1PbW/Pht4dqVP93f3LwDekeU5m24A3Jzux+T6Lu1dwDsr7/ktwGn99XHAtcBNK/mvBO7YX79irW7gRGBvqe90X7b7ubQjgGVgF/Afgc8Ad2j4Hy0D/8rdj7Tr0k4HXu3y/Qh4A/Bp4EYu318Cr8rKfpNucOfvw7/nY/v/6dGu3OeBx/bX3wVOcmkPBL7fX78V+GOXdoP+vR1L9yP+WZcmYO/a96V/trtv++bbNR6247Xdmr1Hmtn/LqRd7K5v0X+Al0j7JJqBy3Nklv8HlTaPoftHt3AL4NGSHuae7QbO69u80syuzto9plLf2+nE8mfRzbofN7NL1xIlPZdOmjiS7styCHVRvNbvD0kaumerdD+C7+j7+G5JhwLvBF5iZstj6rmSTvLZCIfSSU7/3sx+mvXtFEm/754dRPeeW/h/7vqXdAORvrz/v//A1Xkk8E9rCWb2C0mXA0eRfXfMzCT57xKsfwZXNfZxLpilosi7M11MN/PuMbND+9chZvY7ffolpIPm5pV6LwZGtL5j2lzL+w7X5qFmdn0ze13f5mGSrt/YLmb2f4DLgUcA/wEnMvfr2xfQiXaHmdmhwE8ZXX8BXA1cz5VdAm7i0i8GHpz1+2Az+5GZLZvZK83stsC/oVtuPLHQ5QuB36q9pwpX9nX/tSQvhl4MvCbr2/XM7MwNtrPGj+l+GNa4ef9sJK3/n92YTjpIvjvqZof8B/g2dLP4zzbZxx1lLrS8ZnYJ8PfAf5d0SK+Q+U1J9+6zvBd4lqSjJR0GvLBS3VuA50m6a6/JvrWktX/sT4BbubzvBB4m6YG9UuzgXll0tJn9ALgAeKWkg/q17MOYzBnA6+lmpr9zz29IJ4b/C7BL0svoZt5xfAs4WNJDJO2mW6/5NeOpwGvW3pekm0h6RH99H0m37wf8z+jEx9VCOx+lWzcn9J+Df437gcHMzgceTycFnNA//p/AMySd0H/+1+/fx0Zn+DXOBF7av9c9wMvo/n/QLWeeLOlO/dr6tcDnzOz7dDqN35H079TtIT8L+I2s7nvTKVQXirkYvD1PpBOvvkH3q/5+urUcdF+Is4Gv0IlHHyxVYmbvA15D9w/9OfC3dAoi6JQaL+21oM8zs4vpZskX0w2qi4H/yvrn8jjgBOAK4OWkyqcSZ9DNCu8xs2vc87PpviDfohP5fk26FPDv4afAf6b7IfoR3Uzstc9/DpwF/L2knwOf7fsJ3Rfz/XQD9yI6bXLJ+OMM4CRJ13XPjqJTcPlXSZLBzM4BnkynlLqrmV1Ap2B8E93/8Tt0Cq3N8mq6H9MLga/SfQ9e3ffhXOAP6RRTl/T9fWyfdhnwaOB1dFLRccA/ZHWfTL/XvUioX7AHByiSXgtcamZvnHVfZkGv73iCmT1m1n2Zlhi8QbCgzJPYHATBFMTgDYIFZVODV9KDJH2zN0uraYCDINhiNrzm7bcivgU8gE4T+gXgZDP7xtZ1LwiCEpuxsLo78B0z+x6AOp/IR9Bt9Yzluodex2545GSvwNLPidnY7caRMpbZPfjfpzxtjWGtbpdmhee1tJH349OSRFeft50aKT++rtE0V3PheTWtMV9r3Rra2HzKP6Ak3/p1Un4k33D882H6Qa5NVr/maq61a8r/dMcD73N9u/yK0lY5fPHCa842swe11LWVbGbwHkW6T7mX9b3GfUh6Or271Q1+43o85p0PBGBoZYl96L7Eq+7LuTJc2lC+FdfWytBd23q+a1fTMquuzLJLWx6Ofw6wslpoZyXLt7KeNnT5hu65LWefz4r7nq2uXyvLN3BGkHJlBis+T/qd9WVK1wCDa9evl5Zt/PNr0wHm8y1d48sM3fN0gA2uWR8og1+vuOdph/Qrt43+6/Vr89e/Sn0yhr/+NQCfs3Np5bIrVvjMx48qph985D9vxMR102xm8I771RqdZMxOA04DuNltD7fBvp/ZytTiB7ZvZZD++iWD1E9a2U/5wHVr4NIGNv553+/1ql1a0p2sTCmf8v4M1u/9hOHz2SCfwlyNSeVZvwduYCf5Ctdkk3cln9eQlMrkv8lWSiu9n/y+VAawpfVEDVxG9xkwyDq0VscUq0UDhtMU2CE2M3j3ktqIHs26rWkQ7DcYxrKVxeZZsRlt8xeA49Q5rB9EZ4521tZ0KwjmiyFWfM2KDc+8ZrYi6b/Q2ewuAW8zs69vWc+CYE4wYLmqQZwNm/LnNbOP0nmmNLPUfwhL2TpnNREC3AdVWv+SrW3d5UCZAsQtvJI1r8u3a5BWPiyU2TVYL7OaaXqTpZZf/2b9TpZ7iTq2sF4lWwMPfUNpvtKa1QrXeVtWWtfm98W1bN4f/2bdmj7pT/UDKuZLnJ382rZQvrufXtg0YHUOzYi3/ZjNIFh0DGN5P1NYBcEBgRksz9/Y3dnBK0a3V/aR7Jk0iNDAbrd1tOy2jZbyX0mXbzgsiNCZWJSIyk68HjoxbGlkm8aJ1L6dQSrG+z5QEK8t/5zG+8NXt1msJD1W6i5uG1XqronAlojK4993feuqtkRI1inr10tuCzHbKtJavqmWsGK1YNwzS2LmDYIJGLBcscCbFTF4g2ACBjHzSsZudSLsam2LuUWEhkSMTjTHIxWu5yuJwLmG2mufS5rnqoVVwdpqNJ9PGK95HqmkYKGVl1NBZK1ZQflm6tZSbSJwSSvty1u+9eDuU9G/TSud5FtK34Tt+z61D8Zu5p0/79mYeYNgAobqk82MiMEbBBOImbdnn6iZyXteclp1orLXAi+Teud4MXq3e5qfLp4aaTjx2olXw0wkKzkt1MRmf7/kRNvhMBevXZlEBPZibqbJLonUI84D3rmhwUkBMuP/Sr4GB4YRw46GukeNQbxI7TytSppjKBtmZEYZ6uvTuCPoiyjxMpsXYuYNggl05pH5xDF7YvAGwQTMxLLF4A2ChaPbKjrAxWZhLKngmJDsI4zfKtqdRe1IRRl3MkO2WPMODN76alhwOIB8e8itwRm/HQRlR/1RZ3y3XVWytir4kENl/Zt3onmN6a5rjgmF9WvR4T5Pa3HMz+ouOuZD2Rmh5Jjv00oWa2MwxLLN3zw3fz0Kgjkk9yCbB2LwBsEEYubtGZRcqxILp8kiNKRitBehcycD78AwXPXi8NCVyX1zC44JTgzblYnnq+4wORW2jbp8jM1XdR4YlLaK2s66Kjop5HU0ni1VFoHLZVp8e7v7hi2gPF/pPKuaVVYjseYNggWlm3lD2xwEC0dsFdFJUGuOCcNcvkqOu2kQoSERo5d8Wn7ksfP13eWPj/XndGciudc+t2ieAZYKR+SMaJsLonKrM4OVRGgyTbT/SBMtcl5m8vVIWoPmOW/LCtZfZMuKkgPDiMhb0liXfHth9CjYBrpjcEJsDoKFI8TmIFhgQmEVBAtIzLykFla5Kc5AbjvHLaIGiUVUtuYpLI1zn/3k3ClX9zDxEGo9LrZiOeW3h9z1asXCKmnHHwAwyN/E+pdHlTVmyUOoulXU6lVUrLtx/VpcW+drWV+fy5eZ5alhS2nEgX9pIxZW9dhas2JijyS9TdKlkr7mnh0u6RxJ3+7/Hra93QyC2bGmbS69ZkXLz8npQB6+8IXAuWZ2HHBufx8E+yWdM/7GB++4CXBMnhMlfVnS1yV9sqVfE8VmM/uUpGOzx48ATuyv3w6cD7xgYmtaF09zMXWY7g+4MuWjXxNneBt/NhVkUQJtvKi8Oyvjt45Saysn3mdG76s2vsyoeL1+nRogVaIs+G2WZDsos7BaLVVeuGYKx4QGh4ERMb6lD7X+VCImeOd8+bOqtjhiwhY4458OvAk4Y2zt0qHAm4EHmdkPJd20pdKNrnlvZmaXAJjZJa2NBcEisjbzbrj8+AnQ8zjgg2b2wz7/pS31bvsqXNLTJV0g6YKrr7h2coEgmDMMMbTyawv4LeAwSedL+qKkJ7YU2ujM+xNJR/Sz7hFA8ZfCB9c+5naH2D4Lq1wMUX6G6VoFJcurFO+YkPv9JrjqvC9tHpA7cWYwn2+8FhpY16QDK976K8u3WtRe467LvsL1c6bGW1/VxOGiQ0TuC9siXreK5MnRr62BxNNsJQurYtBtWH9PU4y5LtxJdebdI+kCd39a/71vZRdwV+B+wHWBf5T0WTP71qRCG+Es4BTgdf3fD2+wniBYCCbMsJeZ2fGbqH5vX8fVwNWSPgXcEagO3patojOBfwR+W9JeSU+hG7QPkPRt4AH9fRDsl6wZaWzjVtGHgd+VtEvS9YATgIsmFWrRNp9cSLrfdP3rJJV9x9DkRhHu2kcb9xrh5WHW3YIYnZ/057XN3oHBa45XR45+nWykkTsmFOP4Dsta8pLBhvIjS5MjcspHxHpR1zsCJG9v5Kgar9Fl7HU1reSwQC5e+3YKAciyfifi9YgYP14TXYzbCzBY+25MY6ShxLllWvoJ8EQ68Xov8HL604rN7FQzu0jSx4EL6UyM3mJmxW2lNcI8MggaGPGCm4LKBOjz/AnwJ9PUG4M3CCZgloaQnRdi8AbBBNa2iuaNGYQ76S2ssufDFquqwUpSJlkD+3Om8qNfS1s9bs26O7f4SiysCiFSBumW1Iq3sGL8Ohkyp31/7pXLMxKQ24cuSY71mj4sSnWrqPk8KldfxdG/GEGwut1V6OvI2th9Jt5poeSYP+6+ASP9384LMfMGQQPz6FUUgzcIJmCmmHmFsZQ7264nrtPomODFaC9C71ZmYVXw9S2J01B2Wig5LEB6RGxiYZWJwCuFI2JL1lZ5vuYg3KXnjWJqPbj21tU92k5B7B2JmFAQlUtOCrhtpGksrEj/Z/NCzLxB0EAorIJgATFCbE6Ofs0P9Brgj8FxVkcu30g85IImOrfEKml+a1ZQXvtc0jwPMiuoRCudHANbtrBKNM+JRjm3sBpffiQIt9c+J5ZKSQeKZdKjZaaPxlA70bflSBwAb3FYs7Aq+fpWLayWprewwmLmDYKFJNa8QbCghJEGgI/Pm2mdE9EyOd7fWySktSVitCs/zA0u/ImR/tr9Q3INdeq3O17zvFJxTKg5MJQiKwwKojGkRhveiaIex9cn+OuyOFwVZ0tib1VzXOpDodNZ/5pPmSwYaeSi9shpko1ExIQgWEDMRnUi80AM3iCYSIjNQbCQGDHzImDQZGHln9eiBK5f+vXv7sb+rFYi/vk1Z2nbaFfmmODX0yvD8dtBAAO3fkqXfmVnhuJppiPr1wZLrIqTfc1pv8WxoHb0a/lI19y5orBdVel3eoZVYdvI55tmIrVOdJ43YuYNgglEiM8gWFhizYtkHFSwsFpyFlZetF5y+a7Nu1vYRsotsQY+cJn7BfXbQ/k/J0krbBst5RZWTvT2InVuWle08nJ9WKmcYTVwFlHDzAdDJauq2pGujYGym6yqRraX/P7X+L5Vz70ajBeHu7bGb3HVLazW7qcbjP6Y4HkhZt4gmEBsFQXBAhMKK/wxOJVAYwUN80FKj8FJxGgfaCwTyYdJjN/xx+XszjTH/pe2pHlutrBqjKyQOiykZbz2Oj0GthzHtxg9oWZhVXJSqKTVROCNiNqtmuxEjC759laOi23FEMOYeYNgMZnDibcpYsIxks6TdFEfO/TZ/fMIsB0cGBjYUMXXrGiRBVaA55rZbYB7AM+UdFsiwHZwAGGm4mtWtIQ7uQRYi8X7c0kXAUexgQDbwoU7GUn0Z1X552ULK78G9uvffD1dsrhKgmHnnj94T6Lx20ZVCytX967MY6m0zi1ZW0EaadBHGSwZEOXXNYf58to4zdYUhLtqBVWoq7Gd2nldpS2lrbCwMvaDraI+QPCdgc8RAbaDA4VebJ43mlVokm4AfAB4jpn9bIpy+4Jr/+yKlckFgmAescprApLeJulSSdXgYZLuJmlV0qNautQ080raTTdw/8bMPtg/bgqw7YNr3/r217Pdvaibb+d4+WhJ462tRrpr3hKrzWl/kGwPrf+Y5GKR3zoqbRvlQbMTC6vCdlBeh8+36rZ98mNXShZWuXidWlIVrkec2tcvi04KUN/eGVMeKmdYFayoqvmqDgwlZ/xCcO2p2LRi6nTgTcAZxRakJeD1wNmtlbZomwW8FbjIzN7gktYCbEME2A72Z2xzCisz+xRwxYRsv083QY6dBMfRMvPeE3gC8FVJX+6fvZguoPZ7+2DbPwQe3dpoECwc26hVlnQU8G+B+wJ3ay3Xom3+NGXd3FQBtuXOsMrXCgMVLKQSy6t8zey779IsfVvpGVburKsk0FhmYeXFeMaL5zUf4PSI2LbICkosrFKNuQr5Ro91Gn8sbKr1rRzpWgk01hQloWph1ean21p3WeNdeD/gHBWmHIz1te0eSRe4+9P65WIrbwReYGar05yxFRZWQTAJA+pr3svM7PhNtHA88O5+4O4BTpK0YmZ/WysUgzcIGsgPcdnSus1uuXYt6XTgI5MGLsxg8K4ZaSzlYmoiX7kE/4OXW70nYvT6WxmWjtoh9SNO/HwzMco7DySaZ1d+lVYjjSzQWEGkLhlsQOqosOKazeP4riaiqUuoBCcrxfStxfFtMtgY6cP466qRRknshuLRr1t+DM64Tk6BpDPpDJr2SNoLvJzedsjMTt1ovTHzBsEkLDUAnLq42clT5H1Sa94YvEEwEW2rtnmjxOANgha2cc27UXY8uPaahdXoaXzjLayuLa1/IVlsDSsWVp7d2Tp1X/msP37rqGXbCLJohIXIhPl9qwN/MbJgJZqgP88qPduqFv3PW2VRzGcli60NWVi1lak5PZQiBuYWVtqQhRVz6dAbM28QTMJAc+iYEIM3CFo40GfexJ83U995EdSLqQe5H7zVkYVHofvZB+23oXw7fmtnoFwEdmUatnbyOryv73ImF/rjXncNvbWVE8lzC6thycIqjyZYOBZ2K86MavAVrpVpDchdiiBoS5ljwqAgrtcsvnJHhUZyo7R5IGbeIJjEZAurmRCDNwhaOOBnXvmjX1PSo18LR+JkpVYLGubV3FqK8Vppby2Va6GHhSNySkfiACxr/djVpUa/35IDw+hxsT6NYr6SSK2SkwJQjKawEYeBSpkmEZpcc2xjn3d1eA36ZN/epI4pj4DdjJHGdhEzbxC0cMDPvEGwgCi2ioJgcTngxWZhI1ZJLnGdpvVvuo3kLbEOavyRXE0stNJCg8SraPxZV/l606+HV1hf/+ZeRcubtLCqB+H2lk8Fh/dauJNKNMEWT6JWr6JS9MBamdb1dHH9C1McuZgRYnMQLCCb9CraLmLwBkELB/rMK9ZFy9wRwLttJBZNNTnHfaBLiclPmq20VbTbOfPnwb791lFp22gpt8pKgmaXLaxK20NevF6unI/lra9Gg3C7Mt7aquCkABVrpxHxutBQzdG/wcmg+bjYypZSun/m8mQBuac5Iyopd6AP3iBYWGLwBsECEmvejpJjgpeoUpG6zU+3eO4VuVWVd1JwDhCkx8oOnYzmxdxEpM/CpZd8fXMLqyRAd+F61DGhcNbViGPC0KU5Z4akMlI24pvrnxe00Hlaq2NC8QyrvO6lQr6Cb29Xx2C0/QmIGLxBsLjModjcEu7kYEmfl/SVPrj2K/vnt5T0uT649nskHbT93Q2CGdCLzaXXrGiZea8B7mtmv+gDjn1a0seAPwD+zMzeLelU4CnAX9Yq6ow0+qNfs5+yxJkgMcwoiNBkkRV8xAUrOzCUnRRyIw2n3XXicXI8Tqbp9Y4JicFFNbJCm5GG9wFeXnUOEJnm2B8L2+KkAJVoCpV4uM2OCcV2Cs9r7TQ7JlT6sx8ZaUx8K9bxi/52d/8yurgq7++fvx145Lb0MAjmgHmceZt+hyQt9UHGLgXOAb4LXGVma1qevcBRhbL74vNeecUcrvqDYBI24TUjmgavma2a2Z2Ao4G7A7cZl61Q9jQzO97Mjj/s8I3KLEEwW+Zx5p1K22xmV0k6H7gHcKikXf3sezTw45Y68rXceuWFAsX1LyThTlxkwPysq5KjQmnbCGDJBadJLbHGn23V5XNWVW79662tILW4anXGLx0l2+rAUI/KN36dmx8Rq8I2Uvt5VG1lyts++Vp9vJVX1TFhP7KwatE230TSof31dYH7AxcB5wGP6rOdQgTXDvZXjE5XWnrNiBY59gjgPEkXAl8AzjGzjwAvAP5A0neAGwNv3b5uBsHsEL1DfuE1sbz0NkmXSvpaIf3xki7sX5+RdMeWfrUE174QuPOY59+jW/82I2B3/1OV+8+m/rylCrKfOSd+rrq0fMPZ+/p6Bwa/BZRHFiw5Lfgy+XGxXoyuRVYonW9V214qBd7WSMQELzavPy85KeT5msVZ7+hQiSxYtJBKnA9yn9uSH3KarUUkH6l7X1+nPMNqc2Lz6cCbgDMK6f8M3NvMrpT0YOA04IRJlYaFVRC0sLkogZ+SdGwl/TPu9rN0OqSJxOANgklMFo/3SLrA3Z9mZqdtsLWnAB9rybjzERPWPoXKh5GIkpTkLhIxeik5zj+rzz3w1lclLXJXxovKXiPs2yxrhHcn/rxLxXzeh3eXxrfT3bdpm5MjclqcFCATU2sWVg0ODPm3vGD51C6eF8rk+Rp8e5MyUyqdJ2wJXWZmx09X45g2pPvQDd57teSPmTcIWtjmrSJJdwDeAjzYzC5vKRODNwgmsc3+vJJuDnwQeIKZfau1XAzeIJjAZv15JZ0JnEi3Nt4LvJzORwAzOxV4Gd1265t7zf1Kixg+g6Nf15zx0zTvZbRa8pSubBXhrJtyEcd7HyUeRoVto67MeAf80vo3zzew8rp0UFjbplEGszVvyRmfFO9l1OJhlN9XndRLx726LDVn/I14FdWc8Zu2inJn/KUNmuhuQmw2s5MnpD8VeOq09cbMGwSTMJKTTOaFGLxB0MA82jbv+OBdF1oqn0aLtRUwLCTmz1cTy6fJ20aQbh0NElF5vJN+l1YQtUest1xkhUJkwZrDQclJYSRfi5MCZI71tXwNDu81KyhH61ZRrW5vPWWN53BNu0W0r9gcerPGzBsELcTMGwQLyDZvFW2UGVhYdddLWdpq8stWMEwfMd6xsWmruRbYyWEtmueuzHjnAX9E7EhwbWdJ1SzaMvkaYJcTta/1ovbIEbHOkqrgpFBzcS06KZD60ybOCBUrqBbxulnbPBIoe7wPbz0gt0b7MoFuq2j+pt6YeYOggVBYBcEiYqkZwbwwQ21zRlE8rhzf4rTAyxrvcwup9rlF8wyptni3E5W9AUnN+CIpn/3n0yNiG4/B2WQc35qRBiVNdKM4m2qyc9F2/HVVi1wMYlbuT6v2e5pICWnBDZbbRmLmDYJJhJFGECwuseYNggVks44J28XObxW1ZlyjtP7N8g2T86gqFlYN20YAg0KIlOT8Kcu2pJLtofHr3y5faZ27Bc74SRn3vLBtlJepB9curUVrFk3j8236uFja1tOj21AbWPSahdgcBAvL/I3dGLxBMBEDrc7f6N3ZwSuNWDLtwwofTtVJwW0BVSyaDnLibSIqF7aNuntfZrw/72jEhPUtpZK1Vd6/UhDuPCB3UWyu1J2Iyi5PzZ+32be3ZAWVZWvZHmq1sGoWr5MzrLJCB1KUwDX6YGNfkvSR/j7i8wYHDBpa8TUrpvkdejZdmJM1Xk8Xn/c44Eq6U++CYL9kMxETtosmsVnS0cBDgNfQhTgRXXzex/VZ3g68gonBtWsWVl7uaRChIRFlvJg5zMu7crtxwbELmmfIjuVhvBZ4YDWNcEXbXLKqqjgmlNJyEbiIvCfpAAAHxklEQVRFK123sPKV1fI1+PZm903WVpW0Vm1z/bjYNceEdq2zFtxI443A84Eb9vc3pjE+bxDsD8yjwqolSuBDgUvN7Iv+8ZisY9+dD659+eVzuNMdBJOwCa8Z0TLz3hN4uKSTgIOBQ+hm4qb4vH3Yh9MA7nTHgyzX6rqc65ctIjQUjTRWK5+oN+BIROPcKCLRKq83VPLz7dK89nq8CA21WLvjxelp8rU4I7T68+b5mhwLWp0HSvVm+dq1zQ0GJJU+1JlPI42JM6+ZvcjMjjazY4HHAp8ws8cT8XmDAwmz8mtGbHTXCyI+b3Cg0BtplF6zYiojDTM7Hzi/v546Pm8QLCzzJzXvvGPCoGhhVbhpXP96ESJfi/rq0q2ZQpRB0u0c76hfip6Qp9WiCZaOhS2ta/P7jTjtt0ZMqG5cJmvRBoeD/H4jDvPVraLJ51ZZ9n7W0qYdixrOn7I1bJuDYAKy2YrHJTaz5g2CA4dNKKwkvU3SpZK+VkiXpP8h6TuSLpR0l5Yu7fjMW9wqajrDKpebnFVVYdsor3vZW2W5unMf4KKFVcHaKu/qUiVfaauo5KSQ3zefYZU8X78eFZvddZLQamFVsNDK8rVaWKWicqGdvNxGtq5aMfKziafldOBNwBmF9AcDx/WvE+gsFU+YVGnMvEHQgMyKr0mY2aeAKypZHgGcYR2fpbOhOGJSvbHmDYKJGNQVVnskXeDuT+uNk1o5CrjY3a+ZG19SK7TD2mYxKE72Q59xnZIIDYlIturj4WY1+4NXS6LyiIVVcgyOs7CirEUeFDTHVW1zg5NCfr+RfCXf3pG0LT5qtXisTkVDvRGtdNH3uKbJbsWYtLa9rCUYdoVmc2NPzLxB0MA2a5v3Ase4+6K5sSfWvEEwCQNWh+XX5jkLeGKvdb4H8FMzq4rMEDNvEDSwORtmSWcCJ9KtjfcCLwd2A5jZqcBHgZOA7wC/BJ7cUu/ObxUVLay8ENCy/s09hNzzrImSx1F1Xer6sOxq92vmUSd7X1/5FzmxlmrwMBpJq3gsqWFLqeaJRGn9S9u6csQKqlR+CyysWs7EGokSuLZnNu3adxMWVmZ28oR0A545bb0x8wbBJAyYQ5fAGLxBMBGD4fyFCdx5x4ReRzYkFzkLhxCVRGhg4EXdwrYRlJ0WalEWEvHYXS9TJrWQKltYpQ79k50UclrTWqMEFh3wR8TUkgNDyUSLRtG2XKb5rKuGNqtlasTMGwQLTHgVBcECYgarB7jYDDDYJ8dUtpi9mFkSoYFhQRWdBzMrOS3UzrpKNcLjj3sdkGuEx+dboqJ5bnBSGKm7+ayrydZW4+7Xn5fvWzXHzZZYjmZtc6FM6RjYrr4NHWI10+NuSsTMGwQTsa0yxthSYvAGwSQMrLJvPyt23DFhSb0cY7kI7I0QnKyTHEeT1jcoaJhzQ5CS08JSwfABMo1wQYuci7alY2FHRdvpvwilIGTNmufkebmdzR6Jo0webvLhbdU2V/rQdAzsZoiZNwgWEJvoEjgTYvAGQQMW2uYgWEAsFFYJ+9a+ayQBsMevU0Yd+X2ZsgdDyWmhdJ5V3oeStVVeJnFUSOrO1sYF54aSYz6MrsmL+QprvNrauBQKZWS5WVoPVxzrS0vWJkf6kfZrdZSPey2VmYoDXWEVBIuImYXYHASLis2hbbNsBy1HJP0LcDVw2Y41Op490YeZtz/rPtzCzG7SklHSx+n6WuIyM3vQ1nSrnR0dvACSLtjkYV3Rh/2g/XnpwyITZ1gFwYISgzcIFpRZDN5pDqPeLqIPs28f5qMPC8uOr3mDINgaQmwOggVlRwevpAdJ+mYfyvCFO9TmSHhFSYdLOkfSt/u/h21j+8dIOk/SRZK+LunZM+jDwZI+L+krfR9e2T+/paTP9X14j6SDtqsPfXtLkr4k6SOzaH9/Y8cGr6Ql4C/owhneFjhZ0m13oOnTgXwP7oXAuWZ2HHBuf79drADPNbPbAPcAntm/753swzXAfc3sjsCdgAf1J/O/Hvizvg9XAk/Zxj4APBu4yN3vdPv7F2a2Iy/gXwNnu/sXAS/aobaPBb7m7r8JHNFfHwF8cwc/hw8DD5hVH4DrAf9EF//1MmDXuP/PNrR7NN2P1H2Bj9BZKu9Y+/vjayfF5lIYw1lwM+tjwfR/b7oTjUo6Frgz8Lmd7kMvsn4ZuBQ4B/gucJWZrfRZtvv/8Ubg+ax7k9x4h9vf79jJwbuhMIb7C5JuAHwAeI6Z/Wyn2zezVTO7E90MeHfgNuOybUfbkh4KXGpmX/SPd6r9/ZWddEzYUBjDbeInko4ws0v6COSXbmdjknbTDdy/MbMPzqIPa5jZVZLOp1t/HyppVz/7bef/457AwyWdBBwMHEI3E+9U+/slOznzfgE4rtcwHgQ8li604Sw4Czilvz6Fbh26LUgS8FbgIjN7w4z6cBNJh/bX1wXuT6c4Og941Hb3wcxeZGZHm9mxdP/3T5jZ43eq/f2WnVxg04Ux/BbdeuslO9TmmcAldJFK9tJpNG9Mpzz5dv/38G1s/1504uCFwJf710k73Ic7AF/q+/A14GX981sBn6cLLfk+4Do78P84EfjIrNrfn15hYRUEC0pYWAXBghKDNwgWlBi8QbCgxOANggUlBm8QLCgxeINgQYnBGwQLSgzeIFhQ/j8p3JFtYMODFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO8AAADSCAYAAACilJfKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXuwbFV95z/fPveFgkFEEQFBJ+iAGcUMAj6SEDSKSNSq+OBRiMSJM0YrMDrlaypRZ2KiVTOasZgMYUoCZhB0omUYxChFQHQckYeID0JE4wMF8fISROGe7t/8sfe597dW77V7d/c53afP/X2qus7eez27+6xev/V7rCUzIwiCxaM37w4EQTAZMXiDYEGJwRsEC0oM3iBYUGLwBsGCEoM3CBaUGLyApGMl3bZodY9o94WSPtUh36mSPufuTdKvTtDe+ZL+tL7+DUm3uLTvSXrBuHVO0Id3S/pf9fV+km6WtHWt250X62LwSrpK0j1dP2hJh9T/ZJtm0Ld/lPT7Dc/PlHTdWrc/BX8GvG/lpjQozexCM3vhajZsZl8ws6euZp0T9OEnwJXA6+fZj7Vk7oNX0iHAbwAGvHSunWnmAuA1Dc9Pq9PWHZKeBfyKmX153n2ZMxcC/3benVgr5j54qQbGl4HzgdN9gqQ9JP1XSd+XdJ+kL0raA7i6znKvpAckPduLTHXZZHaWdEYtRt0v6buSun6pfwM8T9LBru7DgKcDF41bdz4DenGzvj9R0o2S7pX0JUlPd2lvk/Sjup1bJD2/0MyLgc93eXOSXivpi4W050n6oaTfru//paTLJd1dt/+qQrmmpcIRkm6qv8ePSdrm8v+BpFvrei+R9ASX9hxJ19blrpX0HJf2JEmfrz+Py4F9szavAZ7sv7uNxHoZvBfWrxdJ2s+l/RfgXwPPAfYB3goMgN+s0/c2sz3N7P91aOdO4ETgUcAZwAcl/fqoQmZ2G5X4dVrW58vMbPs0defUZc6jmi0eA/wVcImkrZKeCrwJeJaZ7QW8CPheoap/BdxSSOvalxdR/Tj9npldKemRwOXAR4HHAScDfynpaR2rfBVwPPAkqh++19btHAf8eZ2+P/B94OI6bR/g08CHqD6PDwCflvSYus6PAtdTDdr/TPbjb2bLwK3AM8Z794vBXAevpOcBBwMfN7Prge8Ap9RpPeD3gTPN7Edm1jezL5nZQ5O0ZWafNrPvWMXngc9RietduIB68Nb9OhUnMk9Zt+cPgL8ys2vq93sB8BBwDNAHtgKHS9psZt8zs+8U6tkbuH+C9ld4JXAucIKZfaV+diLwPTP7azNbNrMbgE8Ar+hY54fM7Mdmdjfwf4Aj6uenAueZ2Q31d/sO4Nn1cuolwLfN7G/qNi8C/hH4XUlPBJ4F/LGZPWRmV9f15txP9XlsOOY9854OfM7NYB9l16/nvsA2qgE9NZJeLOnLtWh2L3ACw2JWiU8C+0s6BjgWeATVjLAadXsOBt5Si8z31nUdBDzBzG4FzgLeDdwp6WIvXmbcA+w1QfsrnEX1g/r1rG9HZ307FXh8xzrvcNcPAnvW10+gmm0BMLMHgLuAA/K0mu+7tHvM7OdZWs5ewL0d+7hQzG3w1mvXVwG/JekOSXcA/x54hqRnANuBXwL/oqF4UyjUz6kG1Qo7/6lqLfYnqMTw/cxsb+AyQF36amYPAn9LJS6fBlxsZg9PWPeDpX4CPwTea2Z7u9cj6hkHM/uoma1IKwa8v9DGTcBTury3Aq8EXi7prKxvn8/6tqeZvWGKdgB+TPV+AKjF88cAP8rTap5Yp90OPLrO79N2Uus7fhX42pR9XJfMc+Z9OZUoeDiVCHUEcBjwBeA1ZjagWv99QNITJC3ViqmtwE+p1r5PdvXdCPympCdK+hUq8WuFLVQi50+BZUkvBsY1j1wAvBr4PVIt87h13wicUr+f44Hfcmn/E/h3ko5WxSMlvUTSXpKeKum4+v3/EvgF1efXxGVZvTv7Kmmbey0Vyv8YeD7wR5L+sH52KfAUSadJ2ly/nqVKeTcNHwXOkHRE/d7+DLjGzL5Xv4+nSDpF0iZJr6b6f7nUzL4PXAe8R9KWegn2u1ndR1GJ+k0z8sIzz8F7OvDXZvYDM7tj5QWcDZxa/2r+B+DrwLXA3VQzTa+eCd8L/N9ahDvGzC4HPkY161xP9c8GgJndD/wR8HEqkfIU4JIx+3s1cB/wIzO7doq6z6T6J1sRO3c6UpjZdVTr3rPrum6lVuxQ/UC8j0oiuYNKafTOpgbq9eh9ko7Okr5JNehXXmeUOmlmP6AawG+T9G/q9/lC4CSqwX0H1fcxlROEmV0B/DGV9HI7laR1Up12F9Va+y1UovRbgRPdMusU4Giq/413AR/Jqj8VOGea/q1nFMH4GxNJLwT+0MxePu++zANJj6Mylz3TzH457/6sBTF4g2BBmbe2OQiCCYnBGwQLylSDV9LxtZvcrZLevlqdCoJgNBOveWszwz8BvwPcRqURPtnMvrV63QuCoMQ0IXVHAbea2XcBJF0MvAwoDt5NezzStuy1zxRNrjIT/G4lnhddy3fMp4nqSwsV62ir26Wl78+K+dLnuxJa34Ovr+U9FPMNTTTufmCN+Wyo7urPL/k5D9tDnZx0XvTbj7S77i6Z1OH6mx76rJkd36Wu1WSawXsAldfNCrdR2dwSJL2eOqZy856P5tBXv3m8Vlr+GVSSGlrLNOcb+qfrki8ro0FzPg3KgyApUyif15GWyQZvvzmt9LxKs5HX1f2uhpO6lweN1wAMCmnLffc8Gxj9fmM+W15O87l7e3jHruc7dl3nZVbur7Er6Mr2u5f50t8fUEzf9oR/nsQVdmqmGbxNv1pDw8bMzqVycmeP/Q6ywZgttv2SmzX/cLaX8Te+TMsAa/5RH2rHioM362dx8JZnMF9HacADqNecJjU/r8q4ul2ifw7Qc3WbH+Su8p6y9+rfu3w7u64tL9Nrzqc8X37f9LyXqXVWPqAdnSbdqn/AYBIxbY2ZZvDeRuU0v8KBVJ43QbChMIwdVhab58U02uZrgUPrgOgtVC5t47ocBsFCMMCKr3kx8cxrZsuS3gR8Fliiisn85qr1LAjWCQbsYDAy36yZagM3M7uMKvKjGwIbs8U2S1Zpbdtq/SquZdM1UGmd27ZO9mtJK6xrh+pI1q9ufddZYZX1u6QA82vhTAJMlFku45Cizaf1mj/IfHmvwpq3eE261raOa16fkpTpNZfR8nhr3v46dCNe890Xg2DRMYwdG0xhFQS7BWawY/2N3RkPXsGqmopKNx1F7VazTwc775Co7cXjFjtv0pYXtVvtvCVTUVl0L5bJJUknUidSZm6y9eKoF+NdZ7WUfSbeBlwyKQ2ZinxS2aQ0JEY3lullafX9L7qLzSD63TZdmSkx8wbBCAzYUfApmCcxeINgBAYx8wJYadekUv5J/I87elgVNcqwqh5WuQrWCiJ1SVudp3XXNhc01Ln46fKZk5t7uVa64KWVPM/K9BI1d7NWekj69ZJui7bZCtrmxEMr87Da+f5K3lkNVDPv+ouejZk3CEZgiP46DH2PwRsEI4iZF7BV1janlbvLCbTNXfO1lZ/MSaNZIzwcsVRy0ij3oaxtzjXUhT60xAH4fF5DbblThHPm8O32CqLxUIVtWumS6FsIbEjuxxCbQfR398EbBItI5R45prJmBsTgDYIRmIkd42paZ0AM3iAYQWUq2t3FZoFtGs/2Uwq4H6q6o4dV0dTT4mHVeZ1cMBVZvi6dwMPK21ass4dV4XrIVLTruueC7AdD60W/7vZrUXeZbXZRNAn5IId8DZ6UaQn0L3lYFQIbkvrGMhWJHeNG1MyA9dejIFiH9MPDKggWj5h5YcLAhG6RCcUghaH6XDb3a9q2AV1XUduLj0mZbLlUMiNN4mE15L1VEKm9t1Q+iaTmoebggypptElpSLItmJSSIIfM3tXFK6uprV0VNPcZnBg9podVrHmDYAGpZt7QNgfBwhGmIirxYzUDE4qCT2dtszU+z+voKmqXtMVt2mYrbdXaVndHDysvZw4STXFaxmu8E+1wHvSQpFljQpv3Vkm8Vr6u6OKVlVeSlC/n0dDFaKptcEJsDoKFI8TmIFhg1qPCav31KAjWGSszb+k1CkkHSbpS0s2SvinpzIY8x0q6T9KN9etPRtW77j2s2tevBQ+bjlu/Juvfrvn8OrnlqJFkbZzlK0YcFUxNed0lr6y8vpJX1fA62a8xC30bSvM3vnzLvl6FNW8eiNTFK6tKc+thn1CKSkoLNz9vwIDBdGveZeAtZnaDpL2A6yVd3nCi5hfM7MSulY7skaTzJN0p6Rvu2T6SLpf07frvozu/jSBYMFa0zZPOvGZ2u5ndUF/fD9xMdVDfVHT5OTkfyI8vfDtwhZkdClxR3wfBhqQKxm8dvPtKus69Xl+qS9IhwDOBaxqSny3pa5I+I+lpo/o1Umw2s6vrBj0vA46try8ArgLeNqqurh5W3QPwmzN2PmWhbRvXYj6Xp81zqk28Luxv1RaYkKR5qTD/+S2I1EnA/JDY7E09/nnaiWIdai4/VEfJpDR0YkJzm8OBEr4TPtDB50nLrOQaz1V5ZDD+djM7cmQt0p7AJ4CzzOxnWfINwMFm9oCkE4BPAYe21TepIL+fmd0OlUgAPG7CeoJg3dNh5h2JpM1UA/dCM/vkUBtmPzOzB+rry4DNklrP/V1zbbOk16+IE/0Hfr7WzQXBqmOIgZVfo1A1/X8YuNnMPlDI8/g6H5KOohqbd7XVO6m2+SeS9jez2yXtD9xZyugP19568EFGB21z5+1eS/mGvKXUmJSKxmmhTsEILYdmdw9MKHg+5XV31TYXFK2p51RahoJGOO9DyUurbRvXXr/g5ZV0NC2Taphb6nYfysC9iVJgQ5WvUFkL1XEnUzlpPBc4Dfi6pBvrZ+8EnljVb+cArwDeIGkZ+AVwkln7SJh08F4CnA68r/77dxPWEwQLQZcZtoSZfZERDplmdjZw9jj1jhy8ki6iUk7tK+k24F1Ug/bjkl4H/AB45TiNBsEisbDukWZ2ciHp+eM3Z9jSmE4a7dV1et7pHN/WeN5mBwfrpYU6i9fNStLW7W06bReb1+1jeAsBEPl9u5PG6K1kWx0uClrt/AsqidRtkm4S3FAIbAD3cY8VmCCWBws4eIMggEGcVRQEi4cZ7IiZNwgWjxVT0Xpj5oEJbKoXS6vxYXRc83YyPXUNxm85kqTz2jjZq6qwDmzxnEo6lH+MycHSzW0Ol2m+zs1QyXrWnwbo2hwMrV/ddb/5+ZCHVcH7qpd94F6UTT2xXJ9JGdRP8i1h2zBgOYLxg2AxmTKqaE2IwRsEIzBTzLwI2GkqmsBk1CZqt8bwernXy4VtZZqLJ3JYi0mqTbxO5DpfeeI5Vd4/qmRegm4eW0NlCuJsLlmWPKlaT2NI4nbde3Vt5taunhOHB6lPXFZ3IeihZZztzDbmHlbLg9198AbBghIKqyBYQIwQm0FGb1MuJHWn9dCxVhHYiYxdReiSt1SLV1ZX8TqRExOttNdWZ4UKIvWQ0rSDx5b6WaEkGMFd99NsVoizTTXPWX8K8bidvbIS+bybJ1YhhiNta5yJ1GLmDYKFJNa8QbCghJNGTW/MwIRU19i2laQXjcu1lE4/aI0B7ipqlzTM+WFgvWYNc9spC4moPGguP5Sv5PSRiZ+peN22pU1zs513nEw0zx0dO1qdOfxNs1Z6aFlR3487FuPEhCBYQMygH2JzECwiITYHwUJixMyLBEub+iPztZqEknylhHJ5S/azKqxlSde5xXVy21o2yZfblJpNHEVTE1kwQ5snVmk93PfvO6vbt9VvCXroN/e7q6mIwkfX1Sur/XBtf1PYK6tcuB0bY1+1GRIzbxCMwAiFVRAsKLHmRTI2FTysJhFLSuJ1/rxUt7WZl5K05nxte1gl4mzmQpSYgUruQLkJqGD2GTrErBBna62HaxfE+CGzj6+j2RyjbFVUMvu0HTRW8srK43k9SVdbDvsez7XK1Z/L7OuAmHmDYARhKgqCBSYUVsBSb4rAhLa0rhrmguZ4WNS2Qj6fKS/j8iVn9eZBBs1eWomo3bZdrNc2t2l3S2J4XrfTMPvtYYY+0UTDXNBK51rc0caFhnjeUr5yPG95G5y0zEp88DhLWEMMYuYNgsVkHU68nQ7XPkjSlZJulvRNSWfWz+OA7WD3wCqpqPSaF11kgWXgLWZ2GHAM8EZJhxMHbAe7EWYqvuZFl+NObgdWzuK9X9LNwAFMcMB2T8aWDh5Wntw5qdzP8vo1ra/Z02h4zTv+OtmbE5Iy2XrJm5iSKKDC86H75JoUv5wumaHyIPvkGJMWD6uCSSk5ViU/XqR0GmCLV1Zxr6wsXykyqS0Yf6fH1XgOVotvKpJ0CPBM4BqyA7YlxQHbwcbEhn9M1wOdVWiS9qQ62fssM/vZGOV2Hq69fN+Dk/QxCOaPtbxGUNIbZXkk6UOSbpV0k6RfH1Vvp5lX0maqgXuhmX2yftzpgG1/uPaeT3m8bdm0PLK9rq5oJfE4L1/OV87j70uidm4+6CXmIVc+k/0T0VvN7QwdOF0QqYdmBH+fiMfN3laQBT0UPLSq+kablIa2iy0ECXT1ykq8r1qsZ6XtYocP5G6uq52pFVMreqMbJO0FXC/pcjP7lsvzYuDQ+nU08D/qv0W6aJsFfBi42cw+4JJWDtiGOGA72MjYdAorM7vdzG6or+8HVvRGnpcBH7GKLwN715NikS4z73OB04CvS7qxfvZO4oDtYHdilbTKmd7IcwDwQ3d/W/3s9lJdXbTNX6QsZIx1wLaAzVN4WLWJ021pVshXEo2H05qf5+LnoFBfL/O8T7TSyXavzeL0UP9aRMlE25x4W/k43+y9lupr00qXPKzatpX1ZTrKraUYjqF8iVa6HMywUt/YThftBfaVdJ27P7deLiaM0Bs1fSCtrYaHVRCMwhh2KU3ZbmZHtmUo6I08twEHufsDgR+31bn+HDaDYB1ig/JrFC16I88lwGtqrfMxwH0rptgSM4/n3bo0WtvsyZ3Rk7SO2mZPSVQeEpsL9bWJ2v1Bc748nMxrpb3G2ot7lonaJZF6MBQr3OzMYaUtbIbuW7TSiTOGS0hE5fz0tZK43+IM4su0HEiWaKVpJv//2SlGj7uEnW7NW9IbPRHAzM4BLgNOAG4FHgTOGFVpiM1BMApr2MBgnOLteqOVPAa8cZx6Y/AGwUi0atrm1SQGbxB0YYqZd62Y6eDtydhae1hNsqHXxKYiv04tSC95+dI61x84levxl9z6069zl7K9rpK1sVv79V353HtLiRN+87q26pPrX2I2agl66PtNsbqZlPxJg1Y45BoomoeSp7l1KQky6GhScv/JbWvhnd//2GveMfPPgJh5g2AUlioM1wsxeIOgC7v7zCtgS29MU1FH8Xrggle7mpdaTUUFUdvvwTVsKnKisheHs3wl8To5PDAz05RE6jzO1AdBJJ5c7liEoUkkOf7Pm5fyfLsurfA8/+itsA3r0AHfBXyQQW9IvnbXhTDx0uHa467ahnaQXQfEzBsEoxjtYTUXYvAGQRd295m3J2Pb0hTa5lYV4S65qU1zXPbKSvWSvq1BSduct+NEap9vU5ZvOdFE7yrjRegh/373y5+G3GaeWImY6kXlQbHMQF5DXdYcJ1rppFGXp0WyNX/X5mFVEKkH2QjqqolO65iMaZw01oqYeYOgC7v7zBsEi4jCVBQEi8tuLzYLY+u4pqKWdU1ugtlZpuUs1UHBBDTIYrtK6+RNKpukkrWs+7bzs11L61wV8uT5+s4c1M9mhL7vn7M99ft+/ZsvTP0B1t7bioxdfUo2C2gx+yReWX5tXNhfq05tTsuyDSbwxNqZLTysgmA3YMqoorUiBm8QdGF3n3knCcb3lMRkyEXlfpbWLCr3E9G416mMF69zs9MmJxcu21JjGcjMSO4nfannReMsMMFfe08l5QEM3sOqN/I5pCK1XH8GmVu/f7uDwokHQ1+RD+BPivgghbYDsLvte1XaYna1tooJD6sgWFRi8AbBAhJr3mr/oD16D4/M1+8o7EziLZW0M4HY3E+8rZbSMi7fJmv2toJUY73s2vWHfPUyOa2koV5uCWBYdiJryUMLUpE6FaGzbVOTuF8fwOA/30zUThoa/XzoPtkDK5v+lktitFsiDO1hVecYQ9ssYvAGweKyDsXmLsedbJP0FUlfqw9Jek/9/EmSrqkP1/6YpC1r390gmAO12Fx6zYsuM+9DwHFm9kC9cfQXJX0GeDPwQTO7WNI5wOuoDkcq0tP4ThqTiNBtQQ++vkGb2NxBK705++Z2+Jhid52LwF773Es0z83iNKSidy/RNmfO+oV8XtTe0c/6U9BK9/u5Jtvd9P2ywovT2X9zIl6Xvst8K5/SNNfizFF0FMlPTNg42+CMHBn1wUcP1Leb65cBxwF/Wz+/AHj5mvQwCNYB63Hm7TStSVqqN4u+E7gc+A5wr5mtTKMrhyI1ld15Pu+D9zy0Gn0OgtnSdjbvHGfkToPXzPpmdgTV+SlHAYc1ZSuUPdfMjjSzIx/x6K2T9zQI5sh6nHnH0jab2b2SrgKOoTo/dFM9+448FAkmC0zwtAcpdAtGKK95Mwf/DmvjvD+bCmakzVm+HS7Nr0tLa+Eq364yy96klP1metNRrxD0kK+TSyalnDSeoZtX1qBwZ63zRpv3VXO+rp5Yk06V69HDqou2+bGS9q6v9wBeQHU48JXAK+pscbh2sHExqt+d0mtOdJl59wcukLRENdg/bmaXSvoWcLGkPwW+SnUKWhBsOMT6nHm7HK59E9VJ3vnz71KtfzvTw9jW2zFOkYRW0TiPUS2UK3tblT2svAjtTTa5GcuXadv61Yu6XqQuidMAvURUdiJ0fpJfwUtLBXE6xwc6DInXxTJtYqr/7D2F2OCqQteoP6FiAhG6ZDaMrV+DYDch3CODYAGxmHnpyXhEITChP8k2ngUxOq+rmM+L05noVhK1/TYzOywLTHAi2ma3z0ueryRSe3G6l/U58b4qbDELsOy00g97MdzVnYvaXrze4cVrlcVrX2a5JZ8nPSzNfb5ZpETuE7WTzIuqPJ5KWmh2BjOMOxanMQlJOg84EbjTzH6tIf1YKoXvP9ePPmlm/2lUvTHzBkEXppt5zwfOBj7SkucLZnbiOJXG4A2CUUwZz2tmV0s6ZLW6s8Jq7RISBBuWlXjeNfawenYdufcZSU/rUmAOW7+ONhV1PQqlX1hr5WvcfsGU4fO1rZNLa+OlbG+qJJ9fy1oqc5XWw0uJOSdfTzvvKx+xRJqvl6xZd5V5WLu+am9Oyst42r6Fjof8tUYpuR4kd/5TbZVW/QHoXSORVr6L1Y0q2lfSde7+XDM7d4zabwAOriP3TgA+BRw6qlCIzUEwChv+wcvYbmZHTly92c/c9WWS/lLSvma2va1ciM1B0AFZ+TV13dLjVaviJR1FNS7vGlVu5ntYbdPoPawGHU0PJY+rkjgNmeeUN1e0iNol8bqX/Rp7kdr3LRebSyK1v96k1H7it5L1IvXSkCdW2Utr1/NU1F4uzCq5h1XXNE9pydJOs1fWcIsutasn1spnP66H1XSmoouAY6nE69uAd1HFxWNm51DFCLxB0jLwC+AkMxv5AYfYHARdmGKGNbOTR6SfTWVKGosYvEEwitj6tRK1ugQmtB0U5imJx60nK6hZi5zX1UW8XuqlP8depE7ytWilS2JzrnFP8rV4YrWl7Xyea5szjXVTXcN1dPSqmvpozFJgQzYZJlJmixZ6ghm0MhWtP//ImHmDoAO7vW9zECwkBsp3ql8HzEHb3EFsnlrb3BaY0Hf5vAidlekgXveGzvT1+dQpnxep/XX+3rzYvIR32CgHPeSOHjuftzh2ePIAhq5pnvwgtOnIlza7rm2po5PGUvPjkcTMGwQLyGgnjbkQgzcIOhBr3iBYQOKgMcZY83b02ix57wztR1U0KTWvUSFfvzavjZeyjy81Kfl8mUmpsB72Jpwd5GvZQWNa7r2V3zfR5pWVPm/Z62odiJG25PvXbVtZa7ga3ZCti/ebEzNvEHRh/Y3dGLxBMBID9dff6J2bh1W7F1WzUS2Pud1cLJ3HhzYHFmx21bWJ2j7Ni6y5+OlF6h3uo81NMUkMr/tJ9x5buVdWsg9Ws7Q4TM9fjva8GireqqWZ/e++tXidFbeVzd/rSpENYCrqbISrDxv7qqRL6/s4nzfYbdDAiq95MY4F/UyqY05WeD/V+byHAvdQnc8bBBuStYznnZROso+kA4GXAO8F3lwHDh8HnFJnuQB4NyMO1xbGlhWReAJ/s7btYb2ovDmruxSP68XhIRHYezHh85W9oPrJFjnlExO8lrrksZV7ZZW00r1eWdvc1SsrKe/lw0H271EolscHj0suDvtQVvNLiaVua4QkV6Z9n2SsacGdNP4CeCuwV33/GDqezxsEG4H1qLDqckrgymbR1/vHDVkb350/XPveu9ahd3cQjMJGvOZEl5n3ucBL613ttgGPopqJO53PW++idy7A4U/fYts03vm8iXa4Y74hTbYTo/uF0w+GtsHx+RKRrOzYscN2fZxeizzscOHSEs3zoPE6r69NK10KdEj6kEuf7q3ny4ckW0lL7X+Tp5OgV4XESWNI1K7TxtI2r08njZEzr5m9w8wONLNDgJOAfzCzU4nzeYPdCbPya05ME6/1Nirl1a1Ua+A4nzfYmNROGqXXvBjL0m5mVwFX1ddjn88bBAvL+pOaZ31iAmztYCIaTLDO9fkG2aKulC9ZJ9Oy5i2sjR9uO/0vMRtlAQzeROX6+nDLgjE5CaGjSSn5GFpkrCSYwZfJu9N376NX0F1kX++gN+0eVmW8iSkxL7nPIP9ebcITEzRYf2FF4dscBCOQzVc8LhGDNwi6MEfFVImZx/NumSKquS0wIfFiykSi5HBsvDjrReN8O9SCt5Trfy97LyWROje/lMRjX99S9s+SH9DdCR+Y0OKVVRSv86/Kd6GLCD3ErjJe5M1Fax8Lbe4g8dxas1TY4jcJRsgH3c69rsaJ5wVi5g2CxUS7+8wbBIuJwe6usKq0zaPzlT6mzfl2Mu52sxNN2z7mfkGEHmTn+Cb4AAAGNklEQVSydpo2aHyei7IlkTrP1/OeVEm+Zg+t6t4HSnhPrE3FfEkdbeJwKS3XUJc00U6EHuTLohl5XPlPa6mlzTwIonPlMfMGwWIS2uYgWEQM6O/mYnMQLCbz9WEuMeM9rMSWiQ5bruhnH2CyB1VLuWQZ5wPUvUdUtubteVOG36rVm5DyNXhhPbyUrQP9/la+Ct+fNm+raclPN+xsKiqltXXVBfR7k9DAe0FlpqJN/pQ/bzbKq048rJpnxqFA/zrKaOx/wykUVpLOA1ZCa3+tIV3AfwNOAB4EXmtmN4yqdzUPkgmCjYlRGZlLr9GcDxzfkv5i4ND69XpG7EizQgzeIBiJwaBffo0qbXY1cHdLlpcBH7GKL1PFyu8/qt4Ze1jB1g6nxg0K64vNys05u/KlgQkpXtwuidqDIRG42fSUiMbZb19JpH44y+fT2kxKSZnE+8ptKzsUjF8IMmhjElNRoUw/95Yae4/V7qRicyFIIQvGn9hUtLbB+AcAP3T3K9tK3d5WKBRWQdCF9jXvvpKuc/fn1jvIdKXztlKeGLxBMAoz6LeKx9vN7MgpWrgNOMjdF7eV8szYw0ps7qJFdb9DeWyuJ43NdWLukFZajfkSJ6H8wK6kD+55oq1Ov9CHnUO811Dnv6FJEESibS54R5HFDneV/PxWTqU435YyXbXNyXvtqEUZJIEJ2Va9yX5k5QCGzQUR2Eu4vXzb3Z3a5jHF4LU1FV0CvEnSxcDRwH1m1ioyQ8y8QdABm8pJQ9JFwLFU4vVtwLuo5x4zOwe4jMpMdCuVqeiMLvXG4A2CUVjZjtypuNnJI9INeOO49c48MGHz2Lvr78rfJkInW7kMxfN6p/5m8TrfraUkUpfE6apuv9Wqd+zIYnP9CQwtWukivrpJlLltWmRHrjlOROXS9jZZXf7ANh+00Heicq6RHphL6zVrlPM0L15vdhrmISeN5l6PJtwjg2ABsQgJDIKFxdq1zXMhBm8QjMKmU1itFbM3FU1xolx+2p6n7SBov6YqrY2Ht4v1Hlb+6BIXzNBiXqIQAFGluXZLJqWua9n8bU+wBu538HoDkrXyYNBsHsoPNu9TXtt2we8LNkkAw1KveevXsQMTplBYrRUx8wbBCMwsxOYgWFRsHR40JpthkLGknwI/B7bPrNFm9o0+zL39effhYDN7bJeMkv6eqq8ltptZW8jfmjDTwQsg6bop/UCjDxug/fXSh0Um4nmDYEGJwRsEC8o8Bu84cY5rRfRh/u3D+ujDwjLzNW8QBKtDiM1BsKDMdPBKOl7SLZJulfT2GbV5nqQ7JX3DPdtH0uWSvl3/ffQatn+QpCsl3Szpm5LOnEMftkn6iqSv1X14T/38SZKuqfvwMUlb1qoPdXtLkr4q6dJ5tL/RmNnglbQE/HeqbS4PB06WdPgMmj6f4W033w5cYWaHAlfU92vFMvAWMzsMOAZ4Y/2+Z9mHh4DjzOwZwBHA8ZKOAd4PfLDuwz3A69awDwBnAje7+1m3v7Ews5m8gGcDn3X37wDeMaO2DwG+4e5vAfavr/cHbpnh5/B3wO/Mqw/AI4AbqLZb2Q5savp+1qDdA6l+pI4DLqXywp5Z+xvxNUuxubS95TzYz+o9guq/j5tFo5IOAZ4JXDPrPtQi643AncDlwHeAe81s5WTstf4+/gJ4K7vC9R8z4/Y3HLMcvBNtb7lRkLQn8AngLDP72azbN7O+mR1BNQMeBRzWlG0t2pa0ctTH9f7xrNrfqMwyMGGi7S3XiJ9I2t/Mbq93pr9zLRuTtJlq4F5oZp+cRx9WMLN7JV1Ftf7eW9KmevZby+/jucBLJZ0AbAMeRTUTz6r9DcksZ95rgUNrDeMW4CSqLS/nwSXA6fX16VTr0DWhPkTqw8DNZvaBOfXhsZL2rq/3AF5ApTi6EnjFWvfBzN5hZgea2SFU3/s/mNmps2p/wzLLBTbV9pb/RLXe+o8zavMiqmMjdlDN/q+jWm9dAXy7/rvPGrb/PCpx8Cbgxvp1woz78HTgq3UfvgH8Sf38ycBXqLYc/d/A1hl8H8cCl86r/Y30Cg+rIFhQwsMqCBaUGLxBsKDE4A2CBSUGbxAsKDF4g2BBicEbBAtKDN4gWFBi8AbBgvL/Aa00T21BQB0iAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO8AAADSCAYAAACilJfKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHVZJREFUeJztnX20JGV95z/fvvNy5wUYcUCRGRkiY86AUYwoeHA3hmBENkH/MIlEIywYc85qYtREJPG4viQnms0GzYmRxZUdSVSiiYkjUZEAEyUGBBZlgQkvosAo7zIoY4S53b/9o54eqvt2VT1d1V1dde/vc06de7vqeeuqfur5/p63n8wMx3HaR2fWBXAcpxxeeR2npXjldZyW4pXXcVqKV17HaSleeR2npSyJyitpu6Q/mnCaZ0q6apJpLjckrZH0BUmPSvrsrMuz1GhV5ZW0U9IjklbPuixpqlb0EL8r6bGh4xmTLGdEOY6R9JVwj/dIul7SqRWSfDXwNOCpZvYrEyqmE2hN5ZW0BfhPgAGnzbQw0+HfzGz90PH94UCSVsScK0LS3IjTXwAuI6lwhwK/A/xw3LRT6R8B3GZmC2XScPJpTeUFXg9cDWwHzhhxfaOkyyT9SNK/SDoCQAnnSXogyLcbJT0nXDtI0kWSHpR0l6R3SVp0TyRtkWTpShJUwBskbQPOB14cWss94fpqSX8m6W5J90s6X9KaMl9c0nclnSPpRmCvpBUZ57aFcu2RdLOk01JpbJf0UUlflLQX+PmhPDYCRwIfM7MnwvGvZnZVuL5IXYR7clRG+l8F3g38WrgvZ0t6lqQrJD0s6SFJn5S0IZXeZkmfC8/jYUl/mbp2lqRdQRVc2n++yxoza8UB3AH8N+AFwD7gaalr24EfAf8ZWA18GLgqXHs5cD2wARCwDTgsXLsI+DxwALAFuA04O1w7M5XGFpIWf0Uqz53AG4bDpq5/CNgBHBzS/wLwJxnfbVH8oevfBb4JbAbWjDoHrAz36A+AVcBJ4Z78dOoePQqcSPLSnh/KQ8DtwCXAq9L3N+c7GnBUVvrAe4C/SYU/CnhZeEaHkFTwD4Vrc8C3gPOAdSH+S8K1V4Xvtg1YAbwL+Pqsf5OzPmZegKhCwktChd0YPv878NbU9e3AxanP64Fu+GGfFCrlCUAnFWYOeBw4OnXut4Cdwz/WcStvqAh7gWelzr0Y+E7G9zsTWAD2pI5vp65/FzhrKM7AORKT4r6h7/hp4D2pe3RRwX3eBPwl8G2gFyrX1lHfMZwbrrwXDV0fqLwj8nsVcEPq/jyYvsepcF8ivFTD5w7wY+CIWf82Z3m0RTafAXzFzB4Knz/FYul8T/8fM3sM+AHwDDO7guQH+RHgfkkXSDoQ2EjSQt2VSuMu4PAJlPcQYC1wfZCwe4Avh/NZXG1mG1LHs4au3zMiTvrcM4B7zKyXOjf8fUalsR8z221mbw55H0HyArooL05EGfcj6VBJF0v6nqQfAn9D8hwgedHeZaPt4yOAD6fu5Q9IXpCTeFatpfGVN9iJvwr8nKT7JN0HvBV4nqTnpYJuTsVZTyJXvw9gZn9hZi8AjgGeDfw+8BBJa562nZ4JfG9EMfaGv2tT556e+n94adZDwH8Ax6Qq40Fmtj7mO2cwavlX+tz3gc1DNvvw94leQmZm95C88J4TTu0l9f0lPX1UtIJk/ySEea6ZHQi8jqQSQlLxn5nR+XYP8FtDL7c1Zvb12O+zFGl85SWRVl3gaODYcGwDvkbSidXnVEkvkbQKeD9wjZndI+mFko6XtJLkB/gToGtmXeAzwB9LOiB0gLyNpDUYwMweJKkEr5M0J+ksIN0y3g9sCnkTWr+PAedJOhRA0uGSXj6pmzKCa0i+3zskrZT0UuCXgYtjIkt6iqT3SjpKUid0YJ1F0kkIiT16jKRjJfXt2XE5AHgM2CPpcJKXaJ9vAPcCH5C0TtK8pBPDtfOBcyUdE8p6kKRlP/TUhsp7BvB/zOxuM7uvf5BI4dem3tSfAv47iaR6AfDacP5Akor0CImMfBj4s3Dtt0l+8HcCV4U0Lswox2+S/NgeJmnB02/9K4Cbgfsk9aX9OSSdLFcHifjPwE/nfM9+b3X6eGHejUljZk+QDKG9gqTl/yvg9Wb275FJPEFi2/8zyfDQTSR9AmeG9G8D3heu305yv8blvcDPknRs/RPwuVT5uyQvm6OAu4HdwK+Fa/8AfBC4ONzLm8L3XNYodAA4jtMy2tDyOo4zAq+8jtNSvPI6TkupVHklnSLpVkl3SHrnpArlOE4xpTuslEw8v41kuttu4FrgdDO7ZXLFcxwni7FXo6R4EXCHmd0JIOli4JVAZuVdpXmb17rwKfKl4Z3hyxuNf3HgrEaH+Y/eYzxhP8lNvc/Lf36dPfyDbub16298/FIzOyUmrUlSpfIezuB0uN3A8cOBJL0ReCPAPGs5YWX4jgOz+LKxXonaG5n2YJwWvSUyfpAToa77kPcdUpPE1IkLRyqc0ml3hizDcO3qH18SVUyAh36wwNe/nD0Tc/4Z39mYeXGKVKm8o+7qoidvZhcAFwAcqIPNFvaF2NnmdvqB5T68wUipD6OWqo4gNu10NpOuOMM/ribTG/+lmGmW5b2UM16+i17kti+d0ZP/5hWoMxfSiv8uBvQaKAGrVN7dpOYTk6xIWbR43HHajmHss2zZPCuqvPavBbZKOjLM6X0NyfpVx1ly9LDMY1aUbnnNbEHSm4FLSXTqhWZ288RK5jgNwYB9lOhHmTJVZDNm9kXgi7Hh1enQWbt29MUStl+0/VnCth3MqKJdWjX/ppNjvw5886yOxJz4aZtZeXZq2rZOhVtkc/fz6sU/EwO6DezQrFR5HWc5YBj7lliHleMsC8xgX/Pqbs2VtyO0Zj75P1aK1jWcM4khoGkO+0xzbDdNVXkYOwSTk8+gVI4cUhqQzan4w1I9XFN3nPspuvmzRWaCt7yOU4AB+8wrr+O0DgNveVEHzfdlc87NKNM7W5dUnrB8tRb1ROdK2Cyy5HGObFbsrKx0uIze5uF89kvyvfEmTtLyNm8mnLe8jlOAIboNXPrevBI5TsPot7xZRxGSLlTibuemjOuS9BdhXfyNkn42ply19zbbmggHf7HStIyErUuSTzJ+QxgQoGV6pUstUsiQxjnpDcj7oTD7JflYIwOiW002byfZ7TRrA/tXAFvDcTzwUUas0BvGZbPjFJBMj4xcqTYqvtlXlXi5zOKVJK5ijGSr4A2SDjOze/PS9crrOAWYiX2WW3k3Srou9fmCsBQ2llFr4w8n2YQ+E6+8jlNAMlSUK5sfMrPjKmQRtTZ+mJpt3g62elXGtYj4ZW3HEvGsqp06q67ArHJPc2J95KSqzCGgvLINDAHFpW15dnLfHh6j78MQ+2yqVaXU2njvbXacCLqmzGMC7ABeH3qdTwAeLbJ3wWWz4xRSteWV9GngpSS28W4Sn1orAczsfJJltaeS+Lb6MfBfY9KttfKawOZDlpGytNSLrcRwUGWZDAU7HVZNuyULEyKjZ0poGFpYkJdXhlROyetF+fQ/j3E/I2ze/PhmpxdcN+BN46brLa/jFJC0vOWHiqaFV17HKSBiqGgm1D7Dqrc6yTJaDkdvdTN+cSpL5Qkr2YlI9ymSK3WzyIiSL5uz0hqeLZVxzdJhRvc2j7MgJNkGp3l9u97yOk4BLpsdp8U0cVWRV17HKcBbXhKbrjtffBOi7eFSQ0JjRxlk0ovxm23mZpI7hJMma+VPXvyMVUaL4uTZtiPCAKjbHyrKyX9EEr0G2ryFJRq1FlHSwZIuk3R7+PuU6RbTcWZHv7c565gVMa+T7cCw+8J3Apeb2Vbg8vDZcZYkyWL85lXeQtmcsRbxlSTTvQA+AewEzinMrQPd1RHvi0hJU2popUSUiSumlkrlAWJnUsV4JshJSxlDQHnxMoeQUtfG2zus8mL8qVDW5n1af+K0md0r6dAJlslxGkW/5W0aU++wSjvXXr1mw7Szc5yJY4heA3sWy1be+/vbdEg6DHggK2Daufb6p2y2ycrmuHBVZWrTeqgbR+TMq8we5lzZrMhwGTOshmV7/9o4vc3WzJa3rJDfAZwR/j8D+PxkiuM4zaRnyjxmRWHLm7EW8QPAZySdDdwN/Mo0C+k4s6S1kzRy1iL+wti5CbqrMt5UZXqBWyObK8ZvOpE3qJxsjg2Xltc564FL9DYbYqHXwsrrOA70GvgG9srrOAWYwT5veR2nfSy1oaJSWKc9Nu9En1Xznvtkid63avz48TZveqhIqfOj8xrn+Rqw0MAZVs0rkeM0kJ51Mo8YJJ0i6dbgTGzRWgBJz5R0paQbgrOxU4vSdNnsOAWYqVLLK2kO+AjwMpIN1q+VtMPMbkkFexfwGTP7qKSjSbaD3ZKXbu1bv3YzHCaUkpaTnolVMZ+p5d9w4tf2lohfaqgoJ+0yM6yAhV4lkfoi4A4zuxNA0sUki3vSldeAA8P/BxHhMcFbXseJoKDDqsjR2ChHYsMuPN8DfEXSbwPrgJOLyuSV13EKMAplc5GjsRhHYqcD283sf0p6MfDXkp5jZpkemuqtvIJelmxOEb1Od4pSO82sZXcjKbELbJne5tytbqLW846+NtYztcKWt4gYR2JnEza9MLN/kzQPbCRn0Y/3NjtOAX2bN+uI4Fpgq6QjJa0CXkOyuCfN3YQpx5K2AfPAg3mJumx2nAKqTtIwswVJbwYuBeaAC83sZknvA64zsx3A24GPSXoryfviTLP8tZYz6G2OuAlNnnwxAQm8FHqfq/Yw56aVKa8Hb1yMDJ9EbzNU95hgZl8kGf5Jn3t36v9bgBPHSdNbXscpwAy61YaKpoJXXscpxOc2O04rMbzlTYaKVj75f2ycGEq9GGe9YGEZUHkxQl6cLNs2wuYdd6ioqs/xaeAtr+MU4C4+Hae1uM1bSjbXJYcbNwQ0q99KRXk4LFMzk4ueVRUXLrMMGTOsxp2e1Ivx/FAz3vI6TgE+VOQ4LWbZd1hZ3sKEGH9UZZVL9HSgCeS1P8+K8dtE5W1wsm9W7DY4MT3P6c/jbYMjet7yOk47aWDDG+Vce3PYW2eXpJslvSWcdwfbzvLAwHrKPGZFjBZYAN5uZtuAE4A3hT123MG2s2wwU+YxK2LcndwL9H3x/kjSLpJtPcZ3sN2B3qpiATJVNyZV96OqFn0iZWgEsXZuTJQyi+zzEsydYaXFBSvAWAJDRZK2AM8HrsEdbDvLhSCbm0Z0F5qk9cDfA79rZj8cI94bJV0n6bru3r1lyug4s8dyjhkR1fJKWklScT9pZp8Lp6McbKeda6/etNl6/RynOsOqzAZL40eZqrnThBf9BPapKjfDarTHg0VRshYjpE8vks22OFAhs+2YyiKmt1nAx4FdZvbnqUvuYNtZHlhLO6xItub4DeD/SfpmOPcHuINtZznRxoUJZnYV2SJjPAfbeb3NU3U0VnW2fbXojZDA0yRP2maEGzw/fIMynGNPZIaVRmdZRANnaTRvzpfjNA0Deso+IihyNBbC/KqkW8JkqE8VpenTIx0ngmy/BcXEOBqTtBU4FzjRzB6JGXqtfT2vrRyvt89iJe+M1vBOjapSfxJM086zzA+ZMlzD5cmQxwM9zIu2iy3T27w4nTGJcTT2m8BHzOwRADPL9JTQx2Wz4xRhoF72EcEoR2OHD4V5NvBsSf8q6WpJpxQl6rLZcQpRUctb5CUwxtHYCmAryZTjTcDXgqOxPVmZeuV1nBjyW9giL4ExjsZ2A1eb2T7gO5JuJanM12YlWrPNa9jKCJ0x1YUJ9cy+mmj8mZFzr8qY5CUW4w/assO2ccwODjb647i/g2pdEPsdjQHfI3E09utDYf6R4OZT0kYSGX1nXqLe8jpOEQaqMD0y0tHYpcAvSroF6AK/b2YP56XrlddxYqjY+R/haMyAt4UjitqHiliZcRciZEysz+2ZSGNAdQ3vTFqGT7DY0XN9Ixcm5G78ls4ra1ZV1vDSmPewCSN3w3jL6zhF9GdYNQyvvI4Tw7JveWV0VnbD/5FRpiiVS3VWT1M/NVGbpSkxyyhLRucvYMiIk7MNjmVK6MG0rKxsrjA9clp4y+s4MTTwveqV13EKUMWhomnhlddxIlj2slmCub7NWxAuLr0Sdm6JF+ik7dzahpSmSOyQUEy4vOEgyxs2yrSNU3EWXexv/VrrDKup4C2v4xRh3vI6TntZ7i2vZKzMkM0xcras3Cwlr0vlNJrOEpDJecR6jc9cl5ATP1c2Z4bLjtO/Nq751MRH6C2v48TglddxWojbvIlUWbUiSzbHLEwo9/rrzKCHuQlSOf0dprk5eLRszgjXKyGH88owMMFq0QyrvmyOfz7CK6/jtJfZv4sXEePuZF7SNyR9K+wn+95w/khJ1wTn2n8radX0i+s4M6D6BnRTIablfRw4ycweCw7HrpL0JZJFw+eZ2cWSzgfOBj6al5BkrF65kHmtiLJStIxgrCp765LNZeT9pCV0rGzOCpf3DfLSTn+PXoa8HpbkZWRzYSFnRGHLawmPhY8rw2HAScDfhfOfAF41lRI6TgNoYssbtW+zpLngZOwB4DLg28AeM+s3o6P2oe3H3e+fd+HRH0+izI5TL3m+eWfYIkdVXjPrmtmxJFtWvgjYNipYRtwLzOw4MztuxUFry5fUcWZIE1vesXqbzWyPpJ3ACcAGSStC6ztqH9pFdDDmV9Rv85aJ16n4Sp3V4oOs7xprl5Yh1obuZfQ+5JUty5bNC5c1bJS+VvcMq+AB4cMku0f+bzP7QEa4VwOfBV5oZteNCtMnprf5EEkbwv9rgJOBXcCVwKtDMHeu7SxdjGTT9ayjgJSjsVcARwOnSzp6RLgDgN8BrokpVoxsPgy4UtKNJJtHX2ZmlwDnAG+TdAfwVODjMRk6TtsQYUF+xhHBfkdjZvYE0Hc0Nsz7gT8FfhKTaIxz7RuB5484f2coVDSSMb9i38hrMdK2TbJ5MP8GTs+ZID2L81dXVTZHDxuRHaf/edzfREXZPMrR2PED6UvPBzab2SWSfi8mUZ9h5Tgx5L9/Kzkak9QBzgPOHKdIXnkdp4hieVzV0dgBwHOAnUp60p4O7JB0Wl6nVa2VtyNj7YonosLFplemDGPHmfBgXhMWLVQlelZVxPy2srI5S65HyeYxn2lFyyfX0ZiZPQps3J9XMqLze5V7mx3HodIkjTCc2nc0tgv4TN/RmKTTyhbJZbPjFDGB9bxFjsaGzr80Jk2vvI5TgK/nJbH15ucW9v8fFyfurs3NwJadtO3a9CGl2CGhwTgZw0M5tnA3c9F+Z+hz8fBQls27FFYVecvrOEUYKG+7jxnhlddxImjiAEG9shljXcZQUYxkjJXGZeTwJCRrXUNAczETasegO8FBh/xhn+KhnWGyJPRwelkSejh+P864ZlYTLRpveR0nhuXe8jpOK/GtX2FOPdateLw4XKQsrGuG1aRl6nKaYRUjybN6hIvSyuzJTsvpIUne3b8wIf6ZJkNFzXtm3vI6TgRNfN965XWcIgxU7Jm2dupfmNB5Yv//McyVMDbK9DaXyWcS+cYyifLF0C0xESNNzEKEonwGe4tzwqVkc1pS503S6Icb+1l5y+s4LcQnaThOe3Gb13FaiC9MILEz1s9FDBWVsXNjFzBUXoww2adYtTx10i3hOCZrVlVeWplxhhcmRNjGwzb4vt4cMOZwnZnLZsdpLc2ru155HacQA3WbV3trHirqsTZDNsfIx7KStcwMqUnOgpr0DK2mEbuwIXbm1GCcTipcnNTuDqztHUy7q6UzVBQ9qBecjd0g6ZLw2f3zOssG9SzzmBXjjMi/hWTzrD4fJPHPuxV4hMQ/r+MsSSp6TJgKUbJZ0ibgvwB/TOLiRCT+efvbV34CeA8FzrU7PDnDapgYaRnbC90pIVPLbKOzKI2a5PGke7zLbG+TRZ4Ezlubm6aXSiNvhlXWrKrB80OymRILEyYwSaPI0ZiktwFvABaAB4GzzOyuvDRjn9qHgHfw5L7xTyXSP6/jLAXUtcyjMG6co7EbgOPM7LkkTuv/tCjdGC+BvwQ8YGbXp0+PCDryW6Sdaz/2yGg/RY7TaKo71y50NGZmV5pZ3/v81SReFXKJkc0nAqdJOhWYBw4kaYmj/PMGny0XABz5M+ttXSfpbY5fs1tCAte0MKGMPI+lroUIeZRZpNCLFHP5cji7t3gwXFpep7fByZbd5RYmVO6YKnQ0NsTZwJeKEi2802Z2rpltMrMtJG4arjCz1+L+eZ3lhFn2ERyNpY43DsUeR6m+DjgO+B9FRaoyznsOcLGkPyLR6+6f11maFE/SqOpoDABJJwN/CPycmRXOIx6r8prZTmBn+H9s/7yO01qqdTbnOhqD/f55/xdwipk9EJNovXtY0aNv88bai7G2Xxk7t6rNOk27tAkLFroafyFCrJ3cU5zNm5d2lm2bO8Oqb/OOvfVr+WdtZguS+o7G5oAL+47GgOvMbAeJTF4PfDa4+bzbzHKdkPncZscpQBY3JJRHkaMxMzt53DS98jpODDZ7JTRM7et51yoMFU1TDpcaXqpnSCqWJjgdKzPzKlZql1qMMJR0WiqnZXg3Z51vX2qP9bwNWO6rihynrWi5t7yO004MKnRYTYva1/P2e5uHmeR63lIzrCYggetzNDa9fAZl6vibFcd7UshxLqbIGVYqXowwLOP7MnqsGVaG27yO01aW/U4ajtNKDOguc9nsOO3EXDZ3MNZp9LLAGHuxrK1Xzv1JtYfVhBlSdZE3nDMQLsurX7T9mxMuPTyUYzP3r409FLfcO6wcp5UY4Ps2O04bMeg1z01gzbIZ1nYWCsPFSs4yOy/NjT/XvnKeecxNOL26GPwpZz+vAbGZMkXyOm+zBOriYZ/0AvyMxQxDcfYFGT3Wz8BbXsdpMW7zOk4LMYPuspfNxnyQTqUkb2y4EutQJyGH50o44oqlU+I7laFXcUikG2ny9HJ7pVPyeuDKYNq9DBmebiOHZ3KtDFfH95jgstlxWoj5JA3HaSUGZsu88nYk1kbIv1j5WUZKVpW2nQkI7DKyvjYii9aNlJG9iHWzw1K7lyFpF4dLkSr3oOweLbXHnkTjLa/jtBDzJYGO01psufc2O04rMe+wQsC8Rg/4xNiSZW3FMnZqp6JtPJeztemSYMDGjPthZ9myi+ziDPs1z35O27bpfLLsZI37W1ruHVaO00bMzGWz47QVa+DcZlmNM0ckPQjsBR6qLdPRbPQyzDz/WZfhCDM7JCagpC+TlDWLh8zslMkUK55aKy+ApOsKnDJ5GZZB/k0pQ5tZ4r0qjrN08crrOC1lFpX3ghnkOYyXYfb5QzPK0Fpqt3kdx5kMLpsdp6XUWnklnSLpVkl3SHpnTXleKOkBSTelzh0s6TJJt4e/T5li/pslXSlpl6SbJb1lBmWYl/QNSd8KZXhvOH+kpGtCGf5W0qpplSHkNyfpBkmXzCL/pUZtlVfSHPAR4BXA0cDpko6uIevtwPAY3DuBy81sK3B5+DwtFoC3m9k24ATgTeF711mGx4GTzOx5wLHAKZJOAD4InBfK8Ahw9hTLAPAWYFfqc935Ly3MrJYDeDFwaerzucC5NeW9Bbgp9flW4LDw/2HArTXeh88DL5tVGYC1wP8FjieZILFi1POZQr6bSF5SJwGXkMxgri3/pXjUKZsPB+5Jfd4dzs2Cp5nZvQDh76F1ZCppC/B84Jq6yxAk6zeBB4DLgG8De8ysvxfvtJ/Hh4B38OTagKfWnP+So87KO2oZx7Lp6pa0Hvh74HfN7Id1529mXTM7lqQFfBGwbVSwaeQt6ZeAB8zs+vTpuvJfqtS5MGE3sDn1eRPw/RrzT3O/pMPM7F5Jh5G0RlND0kqSivtJM/vcLMrQx8z2SNpJYn9vkLQitH7TfB4nAqdJOhWYBw4kaYnryn9JUmfLey2wNfQwrgJeA+yoMf80O4Azwv9nkNihU0HJwtGPA7vM7M9nVIZDJG0I/68BTibpOLoSePW0y2Bm55rZJjPbQvLcrzCz19aV/5KlTgMbOBW4jcTe+sOa8vw0cC+wj6T1P5vE3rocuD38PXiK+b+ERA7eCHwzHKfWXIbnAjeEMtwEvDuc/yngG8AdwGeB1TU8j5cCl8wq/6V0+Awrx2kpPsPKcVqKV17HaSleeR2npXjldZyW4pXXcVqKV17HaSleeR2npXjldZyW8v8BZJ2WiLUe9+EAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set into eval mode\n",
    "model.eval()\n",
    "likelihood.eval()\n",
    "\n",
    "# Create 100 test data points\n",
    "# Over the square [0,1]x[0,1]\n",
    "n = 30\n",
    "test_x = torch.zeros(int(pow(n, 2)), 2).cuda()\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        test_x[i * n + j][0] = float(i) / (n-1)\n",
    "        test_x[i * n + j][1] = float(j) / (n-1)\n",
    "\n",
    "with torch.no_grad(), gpytorch.settings.fast_pred_var():\n",
    "    observed_pred = likelihood(model(test_x))\n",
    "    pred_labels = observed_pred.mean.view(n, n)\n",
    "\n",
    "# Calculate the true test values\n",
    "test_y_actual = torch.sin(test_x.data[:, 0]) + (torch.cos(test_x.data[:, 1]) * (2 * math.pi))\n",
    "test_y_actual = test_y_actual.view(n, n) / 8\n",
    "delta_y = torch.abs(pred_labels - test_y_actual)\n",
    "\n",
    "# Define a plotting function\n",
    "def ax_plot(f, ax, y_labels, title):\n",
    "    im = ax.imshow(y_labels)\n",
    "    ax.set_title(title)\n",
    "    f.colorbar(im)\n",
    "\n",
    "# Make a plot of the predicted values\n",
    "f, observed_ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax, pred_labels.cpu(), 'Predicted Values (Likelihood)')\n",
    "# Make a plot of the actual values\n",
    "f, observed_ax2 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax2, test_y_actual.cpu(), 'Actual Values (Likelihood)')\n",
    "# Make a plot of the errors\n",
    "f, observed_ax3 = plt.subplots(1, 1, figsize=(4, 3))\n",
    "ax_plot(f, observed_ax3, delta_y.cpu(), 'Absolute Error Surface')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
